import numpy as np

def main():
    # 读取数据
    data_raw = np.array([x1,x2,x3,y1,y2,y3])
    data_raw = data_raw.T  # 输入原始数据，行数为样本数，列数为特征数

    # 数据标准化
    num = np.size(data_raw,0)  # 样本个数
    mu = np.mean(data_raw,axis=0)  # 按列求均值
    sig = (np.std(data_raw,axis=0))  # 按列求标准差
    data = (data_raw-mu)/sig  # 标准化，按列减去均值除以标准差

    # 提取自变量和因变量数据
    n = 3  # 自变量个数
    m = 3  # 因变量个数
    x0 = data_raw[:,0:n]  # 原始的自变量数据
    y0 = data_raw[:,n:n+m]  # 原始的变量数据
    e0 = data[:,0:n]  # 标准化后的自变量数据
    f0 = data[:,n:n+m]  # 标准化后的因变量数据

    # -----相关矩阵初始化
    chg = np.eye(n)  # w到w*变换矩阵的初始化
    w = np.empty((n, 0))  # 初始化投影轴矩阵
    w_star = np.empty((n, 0))  # w*矩阵初始化
    t = np.empty((num, 0))  # 得分矩阵初始化
    ss = np.empty(0)  # 或者ss=[]，误差平方和初始化
    press = []  # 预测误差平方和初始化
    Q_h2 = np.zeros(n)  # 有效性判断条件值初始化

    # -----求解主成分
    for i in range(n):  # 主成分的总个数小于等于自变量个数
        # -----求解自变量的最大投影w和第一主成分t
        matrix = e0.T @ f0 @ f0.T @ e0  # 构造矩阵E'FF'E
        val, vec = np.linalg.eig(matrix)  # 计算特征值和特征向量
        index = np.argsort(val)[::-1]  # 获取特征值从大到小排序前的索引
        val_sort = val[index]  # 特征值由大到小排序
        vec_sort = vec[:, index]  # 特征向量按照特征值的顺序排列
        w = np.append(w, vec_sort[:, 0][:, np.newaxis], axis=1)  # 储存最大特征向量
        w_star = np.append(w_star, chg @ w[:, i][:, np.newaxis], axis=1)  # 计算 w*的取值
        t = np.append(t, e0 @ w[:, i][:, np.newaxis], axis=1)  # 计算投影
        alpha = e0.T @ t[:, i][:, np.newaxis] / (t[:, i] @ t[:, i])  # 计算自变量和主成分之间的回归系数
        chg = chg @ (np.eye(n) - (w[:, i][:, np.newaxis] @ alpha.T))  # 计算 w 到 w*的变换矩阵
        e1 = e0 - t[:, i][:, np.newaxis] @ alpha.T  # 计算残差矩阵
        e0 = e1  # 更新残差矩阵

        # -----求解误差平方和ss
        beta = np.linalg.pinv(t) @ f0  # 求回归方程的系数，数据标准化，没有常数项
        res = np.array(f0 - t @ beta)  # 求残差
        ss = np.append(ss, np.sum(res ** 2))  # 残差平方和
        # -----求解残差平方和press
        press_i = []  # 初始化误差平方和矩阵
        for j in range(num):
            t_inter = t[:, 0:i + 1]
            f_inter = f0
            t_inter_del = t_inter[j, :]  # 把舍去的第 j 个样本点保存起来,自变量
            f_inter_del = f_inter[j, :]  # 把舍去的第 j 个样本点保存起来，因变量
            t_inter = np.delete(t_inter, j, axis=0)  # 删除自变量第 j 个观测值
            f_inter = np.delete(f_inter, j, axis=0)  # 删除因变量第 j 个观测值
            t_inter = np.append(t_inter, np.ones((num - 1, 1)), axis=1)
            beta1 = np.linalg.pinv(t_inter) @ f_inter  # 求回归分析的系数,这里带有常数项
            res = f_inter_del - t_inter_del[:, np.newaxis].T @ beta1[0:len(beta1) - 1, :] - beta1[len(beta1) - 1,
                                                                                            :]  # 计算残差
            res = np.array(res)
            press_i.append(np.sum(res ** 2))  # 残差平方和，并存储
        press.append(np.sum(press_i))  # 预测误差平方和

        # -----交叉有效性检验，判断主成分是否满足条件
        Q_h2[0] = 1
        if i > 0:
            Q_h2[i] = 1 - press[i] / ss[i - 1]
        if Q_h2[i] < 0.0975:
            print('提出的成分个数 r=', i + 1)
            break

    # -----根据主成分t计算回归方程的系数
    beta_Y_t = np.linalg.pinv(t) @ f0  # 求Y*关于t的回归系数
    beta_Y_X = w_star @ beta_Y_t  # 求Y*关于X*的回归系数
    mu_x = mu[0:n]  # 提取自变量的均值
    mu_y = mu[n:n + m]  # 提取因变量的均值
    sig_x = sig[0:n]  # 提取自变量的标准差
    sig_y = sig[n:n + m]  # 提取因变量的标准差
    ch0 = mu_y - mu_x[:, np.newaxis].T / sig_x[:, np.newaxis].T @ beta_Y_X * sig_y[:, np.newaxis].T  # 算原始数据回归方程的常数项
    beta_target = np.empty((n, 0))  # 回归方程的系数矩阵初始化
    for i in range(m):
        a = beta_Y_X[:, i][:, np.newaxis] / sig_x[:, np.newaxis] * sig_y[i]  # 计算原始数据回归方程的系数
        beta_target = np.append(beta_target, a, axis=1)
    target = np.concatenate([ch0, beta_target], axis=0)  # 回归方程的系数，每一列是一个方程，每一列的第一个数是常数项
    print(target)
    print("Generated Principal Component Regression Equations:")

    for i in range(m):
        # Start with the constant term
        equation = "Y" + str(i + 1) + " = " + str(ch0[0][i])

        # Then add each term for the variable
        for j in range(n):
            equation += " + " + str(beta_target[j][i]) + "*X" + str(j + 1)

        print(equation)


x1 = [105,109,107,105,112,107,107,106,107,113,108,114,117,131]  # 混炼时间
x2 = [103.7,105.1,105.5,105.7,104,104.9,104.7,104.8,105.3,105.1,104,105.2,104.9,103.6]  # 混炼温度
x3 = [226,58,248,225,66,225,47,230,237,246,228,68,261,295]  # 混炼功率
y1 = [59.3,59.3,59.3,59.3,60.7,60.7,62.1,62.1,71.6,71.6,62.5,62.5,58.4,58.4]  # 门尼值
y2 = [64,63,63,63,64,63,63,64,63,63,64,63,63,63]  # 硬度
y3 = [12.65,12.77,12.58,12.48,12.51,12.44,13.28,13.11,12.58,13.12,12.41,12.59,12.54,12.73]  # 硫变

if __name__ == "__main__":
    main()